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ABSTRACT 

High-energy photons from cosmological emitters suffer attenuation due to pair production inter- 
actions with the extragalactic background light (EBL). The collective emission of any high-energy 
emitting cosmological population will exhibit an absorption feature at the highest energies. We cal- 
culate this absorption feature in the collective emission of blazars for various models of the blazar 
gamma-ray luminosity function (GLF) and the EBL. We find that models of the blazar GLF that 
predict higher relative contributions of high-redshift blazars to the blazar collective spectrum result 
in emission that is more susceptible to attenuation by the EBL, and hence result in more prominent 
absorption features, allowing for better differentiation amongst EBL models. We thus demonstrate 
that observations of such an absorption feature will contain information regarding both the blazar 
GLF and the EBL, and we discuss tests for EBL models and the blazar GLF that will become possible 
with upcoming Fermi observations. 

Subject headings: galaxies: active - gamma rays: observations - gamma rays: theory - diffuse radiation 



1. INTRODUCTION 

The Energetic Gamma-ray Experiment Telescope 
(EGRET) aboard the Compton Gamma-ray Observatory 
observed the gamma-ray sky between 1991 and 2000 at 
energies between 30 MeV and - 10 GeV. The EGRET 
gamma-ray sky consisted of 271 resolved gamma-ray 
sources included in the third EGRET Catalog of Point 
Sources (Hartmann ct al. 1999) and the diffuse gamma- 
ray emission comprised of emission from the galaxy and 
from the extragalactic gamma-ray background (EGRB). 
The origins of the EGRB are, as yet, unknown; however, 
since EGRET observed a number of resolved, extragalac- 
tic point sources, it is expected that unresolved sources 
of the same populations comprise sizable contributions 
to the EGRB. 

Of the 271 resolved point sources observed by EGRET, 
93 were identified, either confidently or potentially, as 
blazars (gamma-ray-loud active galactic nuclei) and in its 
first few months of observations, the Fermi Gamma-ray 
Space Telescope has already identified 108 blazars (Abdo 
ct al. 2009). Thus, blazars comprise the largest class 
of identified gamma-ray emitters. As such, unresolved 
blazars are expected to have a significant contribution 
to the EGRB. The exact amount of this contribution 
remains undetermined due to the uncertainty in the dis- 
tribution of blazars in rcdshift and luminosity space, the 
blazar gamma-ray luminosity function (GLF; Padovani 
et al. 1993; Stecker et al. 1993; Salamon & Stecker 1994; 
Chiang et al. 1995; Stecker & Salamon 1996, hereafter 
SS96; Kazanas & Perlman 1997; Chiang & Mukherjee 
1998; Mukherjee & Chiang 1999; Miicke & Pohl 2000; 
Knciske & Mannheim 2005; Dermcr 2007; Giommi et al. 
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2006; Narumoto & Totani 2006, hereafter NT06). Hence, 
to this day, it is still unclear whether the collective unre- 
solved blazar emission comprises the bulk of the EGRB 
or only a small fraction of it. 

In addition to the dependence on the blazar GLF, the 
blazar contribution to the EGRB also depends on the dis- 
tribution of blazar spectral indices at GeV energies. The 
spread in the blazar spectral index distribution (SID) 
determines the fraction of blazars with hard spectra, 
which will contribute most significantly at high energies 
and hence will introduce curvature in the shape of the 
unresolved blazar emission (SS96; Venters & Pavlidou 
2007, hereafter VP07; Pavlidou & Venters 2008, here- 
after PV08). However, as with the blazar GLF, the 
blazar SID is also uncertain due to the low number of 
EGRET blazars. The uncertainty in the blazar SID re- 
sults in an uncertainty in the shape of the collective 
unresolved blazar spectrum. Thus, it remains unclear 
whether blazars can simultaneously account for the high- 
energy emission and the low-energy emission. 

The extragalactic background light (EBL) is composed 
of photons from starlight (at optical and ultraviolet 
wavelengths) and reprocessed starlight (at infrared wave- 
lengths). At observed energies beyond the EGRET en- 
ergy range, photons suffer significant attenuation due to 
pair production interactions with the soft photons of the 
EBL (Salamon & Stecker 1998; Chen, Reyes & Ritz 2004; 
Kneiske et al. 2004; Stecker et al. 2006, 2007; Frances- 
chini et al. 2008; Primack et al. 2008; Gilmore et al. 
2009). Thus, any cosmological population emitting high- 
energy gamma rays will exhibit an absorption feature 
at the highest energies in its collective spectrum. The 
strength of such an absorption feature will depend on 
the distribution of sources with respect to redshift and 
luminosity. If the relative contribution of high-redshift, 
high-energy emitters to the collective emission is signifi- 
cant, the feature will be quite prominent. 

Through the absorption of high-energy photons, inter- 
actions with EBL photons will produce pairs of electrons 
and positrons, which will Inverse Compton scatter EBL 
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photons to high energies. The upscattered photons will, 
in turn, pair produce off of other EBL photons, and the 
process continues until the energies of the resulting pho- 
tons are low enough that pair production is no longer 
efficient. For the collective emission of a cosmological, 
gamma-ray emitting population, the effect of the "elec- 
tromagnetic cascade" emission results in a suppression 
at high energies and an enhancement at lower energies. 
In the case of blazars, which could comprise a sizable 
contribution to the EGRB, predictions of the resulting 
enhancement at lower energies could overproduce the 
EGRB if the high-energy emission is high and/or the 
EBL emission is high (Coppi & Aharonian 1997). Thus, 
in order to fully appreciate the blazar contribution to 
the EGRB, we need to include the effects of high-energy 
attenuation. 

While the effects of the full cascade are outside the 
scope of this paper (see T. M. Venters 2009, in prepa- 
ration), we revisit the absorption feature of the blazar 
cumulative emission. Specifically, we seek to determine 
the possible implications of the observation of the blazar 
absorption feature for our understanding of blazars as 
a cosmological population and the EBL. As previously 
indicated, the absorption feature depends on the GLF, 
the SID, and the EBL model. However, all of these in- 
puts (including the EBL model) remain quite uncertain. 
Thus, we seek to demonstrate that the study of the ab- 
sorption feature could be used to constrain the inputs 
of the collective emission and the information they can 
provide about blazars and the EBL. 

Finally, the study of the blazar absorption feature 
could provide insight into the intrinsic blazar spectrum 
or the possible participation of multiple populations in 
producing the EGRB. In the era of Fermi, improved 
blazar number statistics will provide much stronger con- 
straints on GLF models. If the absorption feature in 
the collective unresolved blazar emission is sensitive to 
the GLF, then future studies of the Fermi observations 
of the EGRB could have implications for blazar spec- 
tra or the relative contributions of multiple populations. 
For instance, if the observed absorption feature is more 
prominent than expected from the favored GLF, then 
one might suspect that blazar spectra break above some 
energy (e.g., spectral cutoffs, which as noted in Abdo et 
al. (2009), have been observed in some cases by Fermi). 
On the other hand, if the observed feature is less promi- 
nent than expected, then one might suspect that other 
gamma-ray sources play a significant role in the produc- 
tion of the EGRB. 

In this paper, we revisit the absorption feature of the 
collective unresolved blazar emission at high energies in 
investigating the effects of the blazar GLF and SID and 
the EBL model. In Section 2, we present the formal- 
ism of the calculation of the collective unresolved blazar 
emission. In Section 3, we discuss the inputs of the cal- 
culation and their uncertainties. In Section 4, we present 
the results of the calculation, and we discuss these results 
in Section 5. 



2. FORMALISM 

We define the blazar GLF per unit comoving density of 
blazars at GeV energies, p 7 (L 7 , z), through the following 



expression: 
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where L 7 is the gamma-ray luminosity in erg s _1 at 
a fiducial energy Ef (or, equivalently, E'j times the 
differential photon luminosity dN-y/dtdE measured at 
Ef), V com is the comoving volume, and pl(ch) is the 
luminosity-independent distribution of blazar spectral 
indices (spectral index distribution, SID). In order to 
determine the blazar luminosity-independent SID, one 
must correct the blazar SID for errors in measurement in 
gamma-ray spectral indices (VP07) and for biases intro- 
duced through determining the SID from a flux-limited 
sample of blazars. The resulting pl{o) is given by 
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where p{a) is the SID corrected for measurement error 
in the spectral indices, 
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is the correction for the sample bias (see Appendix [A] 
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and D is the distance measure for the standard ACDM 
cosmology. For an isotropic distribution of sources, the 
number of objects with luminosities between L 7 and L 7 + 
dL-y and spectral indices between a and a + da residing 
within a spherical shell at redshift z with radial extent 
dz is 
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A blazar of gamma-ray luminosity L 7 at a redshift z 
with a power-law source spectrum defined by the spectral 
index a has a photon flux of (see Appendix iBk 
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where d^{z) is the luminosity distance, defined in the 
concordance cosmology by 

d L (z) = -^(i + z)[ z [n A + n m (i + z') 3 ]- 1/2 dz', (7) 
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and t(Eq, z) is the optical depth due to pair production 
on the EBL for gamma rays of observer-frame energy 
Eq originating at redshift z. The total contribution of 
blazars to the gamma-ray background if we ignore sec- 
ondary emission from cascades of primary gammas due 
to interactions with the EBL can be expressed as 

d?V com 

g?L 7 dz da, 

dzdil 

(8) 

where Ie(Eq) is the intensity of the collective un- 
resolved blazar emission given in units of photons 
s _1 cm _2 sr _1 GeV _1 . 
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3. INPUTS 

3.1. Models for the Extragalactic Background Light 

The EBL intensity at the present epoch (z — 0) pro- 
vides an integral constraint on the history of electromag- 
netic energy release in the universe since recombination. 
Measurement of this cumulative output however, cannot 
address its evolution and thus cannot be related to is- 
sues such as the history of star and element formation. 
For this reason, several models have been developed to 
calculate the EBL luminosity density as a function of 
rcdshift, C (v, z), from fundamental astrophysical princi- 
ples. These models were composed with varying degrees 
of complexity, observational constraints, and data inputs. 

The calculation of the high-energy absorption feature 
in the blazar component of the EGRB requires a model of 
the EBL and its evolution over cosmic time. The great 
degree of uncertainty associated with the EBL models 
and their predictions renders selecting the "best" candi- 
date model impossible. Therefore, we use several models 
with widely ranging predictions in order to bracket the 
range of possible EBL realizations. 

Kneiske et al. (2004) treat the EBL-modeling prob- 
lem using separate approaches at the UV-optical and 
infrared wavelengths. In determining the EBL at var- 
ious wavelengths, they make use of a cosmic chemical 
evolution model at the UV-optical wavelengths and a 
backwards evolution model for the infrared (see Hauscr 
& Dwek (2001) for a complete review of the different 
types of models). Additionally, this hybrid model was 
parameterized in terms of the main observational uncer- 
tainties such as the redshift dependence of the cosmic 
star formation rate and the fraction of UV radiation re- 
leased from star forming regions. Thus, the Kneiske et 
al. EBL model consists of several flavors that allow for 
the inclusion of various EBL scenarios that are consistent 
with the available data. Specifically, the Best-Fit model 
best interpolates the data with the important caveat that 
the assumed complete UV absorption by interstellar gas 
introduces a sharp cutoff at 0.1/im. In the Stellar-UV 
model, all the UV radiation produced by the stellar pop- 
ulations escapes to the intergalactic medium after repro- 
cessing by the interstellar gas, with the High- Stellar-UV 
model allowing for a strong UV-field at high redshifts. 
Since the 7-ray sources likely responsible for the EGRB 
at GeV energies are particularly sensitive to the EBL 
density at UV wavelengths, for the purposes of this anal- 
ysis, the Best-Fit and High-Stellar- UV models are used 
to bracket the possible ranges of attenuation. 

Primack et al. (2008) have pioneered the use of semian- 
alytical models that attempt to reproduce the process of 
structure formation and evolution through simulations. 
Recent iterations of this model incorporate highly precise 
knowledge of the local luminosity density at optical-UV 
(Gilmore et al. 2009) and NIR (Primack et al. 2008) 
wavelengths and a well-established cosmological model. 
The key parameters in their approach (those that gov- 
ern the rate of star formation, supernova feedback, and 
metallicity) have been adjusted to fit the local galaxy 
data. With respect to estimates by Primack et al. from 
previous years, this version of the model yields a lower 
luminosity density at optical wavelengths, thereby result- 
ing in a reduced EBL density. Recent TeV observations 
of nearby blazars seem to support such low values (Aha- 



ronian et al. 2006). 

Finally, we consider the most recent EBL model by 
Stecker et al. (2006). In this model, Stecker et al. cal- 
culate the EBL at infrared and optical-UV wavelengths 
separately. At infrared wavelengths, they use a back- 
wards evolution model based on observational knowledge 
of: (1) luminosity-dependent galaxy SEDs, (2) galaxy lu- 
minosity functions, and (3) parameterized functions for 
luminosity evolution. For optical-UV wavelengths, they 
consider the rcdshift evolution of stellar populations with 
an analytical approximation to the more sophisticated 
SEDs used in Salamon & Stecker (1998). The SEDs 
adapted from Bruzual & Chariot (1993) reflect stellar 
population synthesis models for galaxy evolution and the 
observational fact that star forming galaxies are "bluer" 
(brighter in the blue part of the optical spectrum) at 
z > 0.7. Notably, the UV spectra for all SEDs are as- 
sumed to cut off at the Lyman limit, and the effects of 
extinction by dust are not included in the model. The 
former is a matter of debate since it is not really known 
how much UV radiation short of the Lyman limit can 
escape from the star-forming regions, while the latter 
would result inexorably in an overprediction of the UV 
photon density and, consequently, the optical depth at 
higher redshifts. In a similar vein, Franceschini et al. 
(2008) also employ a backwards evolution model based 
rather detailed observations. However, their determina- 
tion of the EBL departs significantly from that of Stecker 
et al., particularly at the optical and UV wavelengths to 
which GeV photons are most sensitive. The differences 
in these models are likely due to differences in the treat- 
ment of galaxy evolution. 

EBL attenuation is a function of the observed 7-ray en- 
ergy E n and the rcdshift z of the emitting source. The at- 
tenuation is generally parameterized by the optical depth 
r (Eo, z), which is defined as the number of c-fold reduc- 
tions of the observed flux, F ^ s , as compared with the 
emitted source flux, i^mittod, at redshift z: 

Eohs ^ ^ ' ^ -^emitted • (9) 

The optical depth is calculated from physical princi- 
ples. Using the cross section for pair production a, and 
assuming isotropic background radiation with spectral 
density n (e) at energy e, the absorption probability of 
7-rays per unit path is given by 
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where 9 is the scattering angle for the 7 — 7 collision, 
e t h = 2to 2 c 4 / [E (1 — cos6>)] is the energy threshold for 
the reaction, and m is the electron mass. Since blazars 
are the sources being considered, redshift is a good choice 
to measure the distance, with the total distance being the 
look-back time (times the speed of light, c) 
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where Hq, CIm, 
parameters. 



and J7a are the well known cosmological 
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Using the expressions above, the optical depth can be 
written as a function of the observed energy Eq and the 
rcdshift of the emitting source 



„ = d _,dl dT(E'.z') 
dl In dz' dl 



(12) 
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dz'— / smO'dO' / den(eV) a(E',e',6') 

where the primed variables (E',e' n(e',z'), 9') refer to 
the values calculated in the comoving frame at z = z'. 

3.2. The Blazar Gamma-ray Luminosity Function 

While the blazar GLF is probably one of the most 
studied and debated properties of the blazar popula- 
tion, to this day, it remains uncertain. Model GLFs are 
typically constructed from luminosity functions in other 
wavelengths (most notably, radio or X-ray) exploiting a 
possible association of gamma-ray emission with emission 
in these wavelengths while correcting for possible differ- 
ences in sizes of emission regions between wavelengths. 
However, since blazar gamma-ray emission is also, as yet, 
not well understood, it is unclear which lower energy lu- 
minosity function(s) could be applicable to gamma-ray- 
loud blazars. Typically, a model for blazar emission is 
assumed and, hence, a functional form of the luminos- 
ity function is adopted. For instance, if one assumes the 
synchrotron self-Compton model for blazar gamma-ray 
emission, then one would expect that the low-energy syn- 
chrotron emission in blazars would be closely related to 
the gamma-ray emission. The unknown parameters (e.g., 
normalization due to relativistic beaming, the faint-end 
slope of the luminosity function) are subsequently fitted 
to gamma-ray data (see e.g., SS96; NT06; Giommi et al. 
2006). 

While such a procedure represents the best that can be 
done with current data, it should be noted that a great 
degree of uncertainty remains as several issues remain 
unresolved. Since it is more difficult to observe fainter 
objects (which are crucial in the calculation of the col- 
lective unresolved blazar calculation), any flux-limited 
sample will be biased toward brighter objects. Thus, un- 
certainties will always be larger in the faint-end slope 
of the luminosity function. Additionally, there is some 
indication that not all blazars are necessarily explained 
by the same emission process and that different types of 
blazars (i.e., BL Laccrtac-likc objects and flat-spectrum 
radio quasars) could form separate populations with re- 
spect to emission (Sikora et al. 2002; Bottcher 2007) 
and, hence, require separate luminosity functions (Miicke 
& Pohl 2000; Dermer 2007). There is also the possi- 
bility that flaring blazars (and different types of flaring 
blazars) could also form separate emission populations 
with respect to emission and require separate luminosity 
functions (SS96) . Finally, due to the large positional er- 
ror circles, there are many unidentified EGRET sources 
which could also have a sizable contribution to the EGRB 
(Pavlidou et al. 2008). A number of these unidentified 
sources could, in fact, be (and are fairly likely to be) 
unidentified blazars. Thus, one might underestimate the 
normalization of gamma-ray blazars with respect to low- 
energy blazars in not accounting for the possibility of 
not being able to identify some resolved blazars; in ad- 
dition, the existence of resolved but unidentified blazars 



would introduce uncertainties to the redshift distribution 
of resolved blazars, which is one of the constraints that 
luminosity functions are typically required to satisfy. 

With the availability of Fermi data, many more blazars 
will be observed, and at least some of the aforementioned 
uncertainty will be alleviated. However, for now, with 
these caveats in mind, for the purposes of this paper, we 
perform our calculations using the best-guess pure lu- 
minosity evolution (PLE) and the luminosity-dependent 
density evolution (LDDE) models of NT06. It should be 
noted that the functional form of the PLE model of the 
blazar GLF originated from radio data, while the func- 
tional form of the LDDE model originated from X-ray 
data. Thus, if the calculated blazar absorption feature is 
sensitive to the blazar GLF, then the observation of such 
a feature could be used as an additional constraint to the 
preferred luminosity functional form and, by extension, 
to blazar emission models. 

3.3. The Spectral Index Distribution 

The unresolved blazar contribution to the EGRB is not 
just a question of magnitude, but also of spectral shape, 
and the spectral shape is sensitive to the distribution of 
blazar spectral indices at GeV energies. If all blazars had 
the same spectral index then the spectrum of unresolved 
emission would be simply a power law. If, on the other 
hand, the SID has some spread, the spectrum will have 
some curvature (SS96; PV08). If the spread is small, 
then even if the blazar contribution dominates the EGRB 
at lower energies, it may not be enough to explain the 
emission at higher energies (PV08). Thus, in order to 
answer the question of the unresolved blazar contribution 
to the EGRB, the blazar SID has to be determined and 
the resulting shape calculated. 

Obtaining the SID of blazars is complicated by the 
presence of large errors in measurement of individual 
blazar spectral indices. If these errors arc not properly 
taken into account, sampling of the underlying intrinsic 
spectral index distribution (ISID) will be contaminated 
by exaggerating its spread, leading, in turn, to exagger- 
ation in the curvature of the unresolved collective emis- 
sion (VP07). Furthermore, the presence of spectrally 
distinct populations of blazars (e.g,. flaring versus quies- 
cent, BL Lac objects versus flat-spectrum radio quasars) 
can also contaminate the ISID. In order to determine 
the ISID of the collective blazar population while min- 
imizing the contamination due to measurement errors, 
VP07 performed a likelihood analysis fitting the third 
EGRET data set of confident blazars to a Gaussian ISID. 
They determined that the maximum-likelihood Gaussian 
ISID can be characterized by a mean (a ) of 2.27 and 
a spread (<To) of 0.2. They also performed the analy- 
sis by dividing the sample of confident blazars into their 
subpopulations, flaring versus quiescent and BL Lac ob- 
jects versus FSRQs. In the case of flaring and quies- 
cent blazars, they found no evidence that the subpop- 
ulations are spectrally distinct (though the lack of ade- 
quate time resolution made dividing the subpopulations 
difficult). In the case of BL Lac objects and FSRQs, 
they found a marginal lcr separation between BL Lac 
objects (a = 2.15, er = 0.28) and FSRQs (a = 2.3, 
(T = 0.19). They also found that the flaring and quies- 
cent blazar populations are spectrally consistent. 

In PV08, the shapes of the unresolved emission were 
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calculated for the collective blazar population and BL 
Lac objects and FSRQs. In the cases of the collective 
population and FSRQs, the curvatures of the shapes were 
not enough to allow the populations to explain all of the 
EGRB, though in the case of BL Lac objects, the cur- 
vature was enough to, in principle, allow BL Lac objects 
to explain the EGRB. However, in all cases, the normal- 
izations of the emission were not determined, and the 
uncertainties in the shapes resulting from the uncertain- 
ties in the likelihood analysis are considerable. 

It should be noted that since the PLE and LDDE GLF 
models do not distinguish between the subpopulations of 
blazars, for the purposes of self-consistency, we also do 
not distinguish between them with regards to their SIDs. 
Thus, for the purposes of this analysis, we include only 
the collective blazar population ISID of VP07(correcting 
for biases introduced in sampling a flux-limited catalog). 
Notably, Fermi has already provided evidence that BL 
Lac objects and FSRQs are spectrally distinct (Abdo et 
al. 2009) . However, since the Fermi blazar catalog is not 
yet complete, it is currently premature to construct lu- 
minosity functions (especially those that distinguish be- 
tween BL Lac objects and FSRQs) from Fermi data. In 
recognizing the importance of correctly treating spectral 
distinctions among subpopulations of blazars, we will re- 
turn to this issue in a future publication. 

4. RESULTS 

The blazar contributions to the EGRB as calculated 
assuming two separate models of the blazar GLF and 
several models of the EBL are plotted in Figure 1 . The 
black lines represent contributions determined assuming 
the LDDE model of the blazar GLF, while the gray lines 
represent contributions determined assuming the PLE 
model of the blazar GLF. For comparison, the blazar con- 
tributions assuming no absorption (solid lines) and the 
Strong et al. (2004) determination of the EGRET EGRB 
(data points with statistical error bars) are also plotted. 
Since the GLFs used include the maximum-likelihood pa- 
rameters determined by NT06, the blazar contributions 
comprise ~ 50% of the overall background 6 . As demon- 
strated in SS96 and PV08, there is considerable curva- 
ture in the unresolved blazar emission due to the spread 
in the blazar SID indicating the increasing importance of 
blazars with harder spectral indices at higher energies 7 . 

The most striking feature in Figure 1 is that of the sup- 
pression at high energies due to the considerable amount 
of absorption by pair production interactions with EBL 
photons. Certain EBL models are quite distinguishable 
from the others. Most notably, the Kneiske et al. (2004) 
high UV model and the Stecker et al. (2006) model pre- 

6 However, as discussed in Section 3.2, the determinations of 
the GLFs and their parameters remain highly uncertain. In the 
likelihood analysis performed by NT06, the most likely level of un- 
resolved emission from blazars is ~ 25% — 50% of the EGRB. Nev- 
ertheless, in several of the cases presented in NT06, parameters for 
which unresolved blazars can account for 100% of the background 
are within the 1<t contours. 

7 There is some uncertainty in the determination of the param- 
eters of the blazar spectral index, which will result in uncertainty 
in the overall shape of the spectrum of the unresolved blazar emis- 
sion. As indicated in PV08, this uncertainty in the spectral shape 
is quite large for EGRET blazars, but will improve considerably 
with Fermi observations. For this reason, and because Fermi is 
currently taking data, we simply calculate the EBL absorption for 
the best guess spectra. 
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Fig. 1. — Collective gamma-ray emission of unresolved blazars 
as a function of observed energy calculated for PLE and LDDE 
models of the blazar GLFs and several models of the EBL. Black: 
the collective emission for blazars assuming the LDDE GLF. Gray: 
the collective emission assuming the PLE GLF. Solid: the collec- 
tive emission assuming no absorption. Thin dashed: the collective 
emission including absorption assuming the Kneiske et al. (2004) 
best-fit EBL model. Heavy dashed: the collective emission includ- 
ing absorption assuming the Kneiske et al. (2004) EBL model with 
the high UV component. Dotted: the collective emission including 
absorption assuming the Gilmore et al. (2009) EBL model. Dot 
dashed: the collective emission including absorption assuming the 
Stecker et al. (2006) EBL model. Double dot-dashed: the collec- 
tive emission including absorption assuming the Franccschini et al. 
(2008) EBL model. 



diet a greater degree of absorption than the other three 
models. This is due to the fact that the Kneiske high UV 
model and the Stecker model predict a higher amount 
of UV background radiation than the others. Since the 
pair production cross section as a function of the center- 
of-mass energy peaks at the electron mass, one would 
expect that gamma-ray photons of energy ~ tens of GeV 
are most likely to interact with UV background photons. 
Thus, unsurprisingly, models with high UV backgrounds 
will result in more suppression at high energies. 

Another striking observation from Figure 1 is that the 
high-energy suppressions for the PLE model are consis- 
tently steeper than those of the LDDE model. The differ- 
ent appearances of the features can be explained by con- 
sidering that the blazar GLF is the distribution of blazars 
in luminosity and redshift space. Since the PLE model 
suppressions are steeper than those of the LDDE model, 
one would conclude that high-redshift blazars contribute 
more to the high-energy emission in the PLE model than 
in the LDDE model. In investigating this possibility, we 
plot the unresolved emission evaluated at several energies 
as a function of redshift (Figure 2). As in Figure 1, the 
black lines represent the emission assuming the LDDE 
model of the blazar GLF, and the gray lines represent 
the emission assuming the PLE model of the blazar GLF. 
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Fig. 2. — Collective gamma-ray emission of unresolved blazars 
as a function of redshift calculated for the PLE and LDDE models 
of the blazar GLF and evaluated at several energies (not including 
absorption). Black: the collective emission for blazars assuming 
the LDDE GLF. Gray: the collective emission assuming the PLE 
GLF. Solid: the collective emission as a function of z evaluated at 
an observed energy of 100 McV. Dotted: the collective emission as 
a function of z evaluated at an observed energy of 1 GeV. Dashed: 
the collective emission as a function of z evaluated at an observed 
energy of 10 GeV. Dot dashed: the collective emission as a function 
of z evaluated at an observed energy of 100 GeV. 

As can be seen in Figure 2, the emission for the LDDE 
model peaks at lower redshifts (z <~ 0.05) than that for 
the PLE model (z ~ 0.6) 8 . Thus, high-energy photons 
in the PLE model suffer more attenuation due to inter- 
actions with the EBL photons than in the LDDE model. 
Additionally, the emission is more sharply peaked in the 
LDDE model than in the PLE model indicating that the 
participating blazars in the LDDE model are more con- 
centrated to a particular epoch while those in the PLE 
model are more spread out over epoch. 

5. DISCUSSION AND CONCLUSIONS 

We have calculated the blazar contribution to the 
EGRB including attenuation of high-energy photons due 
to interactions with the EBL. We have found that (a) the 
EBL attenuation of high-energy gamma rays results in a 
suppression in the spectrum of the unresolved emission; 
(b) the shape of the high-energy suppression depends on 
the blazar GLF model and the EBL model; and (c) the 
high-energy suppression for the PLE model is steeper 
than that for the LDDE model indicating that highcr- 
rcdshift blazars contribute more to the emission at high 
energies in the PLE model than in the LDDE model. 

As demonstrated in this paper, the suppression at high 
energies can be a probe for the underlying redshift 
luminosity distribution of gamma-ray blazars. If, in fact, 
blazars do comprise the bulk of the EGRB, then this 

8 In both cases, the redshift where the emission peaks increases 
as the minimum luminosity included in the integration increases. 
This is to be expected since high-luminosity objects that contribute 
to the EGRB will be distributed toward higher redshifts. The 
inclusion of low-luminosity blazars allows more low-rcdshift objects 
to participate. 



suppression should be observed in the EGRB that will 
be measured by Fermi. Thus, the observation of the EBL 
absorption feature at high energies could provide infor- 
mation about the relative contribution of blazars to the 
EGRB. 

However, it is conceivable that any features observed 
in the EGRB could also be attributable to other ef- 
fects. Multi-wavelength observations of blazars reveal 
the characteristic doubly-peaked synchrotron/ Compton 
structures in blazar spectra. Thus, while it is true that 
when viewed at small energy intervals in the GeV band 
blazar spectra appear as single power laws, we do not ex- 
pect this observation to remain true as the high-energy 
end of the observational energy interval increases. In- 
stead, blazar spectra are expected (and, as noted in Abdo 
et al. (2009), in certain cases have already been observed 
by Fermi) to break beyond GeV energies. Such spectral 
breaks will also manifest as a feature in the collective 
emission of blazars. In this case, the feature would be 
sensitive to the nature of blazar spectra at high ener- 
gies, which, in turn, encode information about blazar 
emission mechanisms. There is also the possibility that 
blazars could account for the bulk of the emission at 
lower energies, but not at higher energies. If the collec- 
tive gamma-ray emission of another gamma-ray source 
peaks at high energies, then we would also expect to ob- 
serve a high-energy feature due to the transition to a 
different population. An intriguing example of this sce- 
nario would include high-energy emission from dark mat- 
ter annihilation, which peaks near the mass of the dark 
matter particle (e.g., Ando et al. 2007; Siegal-Gaskins 
& Pavlidou 2009). In order to distinguish between the 
EBL suppression feature and features arising from other 
possible scenarios, it will be essential to constrain (to the 
extent possible) the expected shape and strength of the 
high-energy absorption feature in the blazar EGRB con- 
tribution, as well as to quantify the associated theoretical 
uncertainties. 

Due to its substantially increased sensitivity with re- 
spect to that of EGRET, Fermi will resolve more than an 
order of magnitude more blazars than EGRET did. The 
number of blazars that Fermi will observe will also pro- 
vide insight into the blazar GLF (see, e.g., SS96, NT06). 
The measured GLF will allow us to not only deter- 
mine the shape of the blazar absorption feature, but also 
the collective blazar emission as a function of redshift. 
With the -Fermi-measured blazar GLF, we can compare 
the rcdshift-dcpcndcnt gamma-ray emission with that of 
blazar GLFs determined from lower-energy observations. 
In so doing, we can determine how closely the gamma- 
ray emission in blazars is tied with the emission at lower 
energies. 

Fermi measurements of blazar spectral indices would 
also allow further investigation into the possible existence 
of a spectral subpopulation of blazars with harder spec- 
tral indices such as high-frequency-peaked BL Lac ob- 
jects (HBLs). There is some speculation that the harder 
spectral indices of the HBLs could be indicative of a 
greater contribution due to a particular emission mech- 
anism relative to the contribution due to the emission 
mechanism that is more prevalent for blazars with softer 
spectral indices 9 . So far, HBLs have mainly been ob- 

9 For instance, in the leptonic emission scenarios, the harder 
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served at low redshifts, and as such, one would expect 
the absorption feature of the collective spectrum to be 
not very prominent. However, if Fermi measurements 
reveal that there are more hard blazars at high redshifts 
relative to soft blazars than indicated by EGRET data, 
then the effect of their intrinsic spectra would be to flat- 
ten the overall collective spectra of blazars. Moreover, 
the collective spectrum of their subpopulation would ex- 
hibit quite a prominent absorption feature, and hence, 
their high-energy emission would further impact the col- 
lective blazar spectrum through electromagnetic cascade 
radiation (T. M. Venters 2009, in preparation). Further- 
more, investigating the absorption features of distinct 
spectral subpopulations of blazars could provide intrigu- 
ing insight into blazar gamma-ray emission and their evo- 
lution with cosmic time. 

Thus, the study of the collective unresolved blazar 
emission from Fermi observations should provide impor- 
tant constraints about both blazars and other possible 
sources of gamma-ray emission. 

In this paper, we computed only the attenuation of the 
high-energy emission due to interactions with the EBL. 
However, when high-energy photons interact with EBL 
photons, they initiate electromagnetic cascades, which 
generate emission at lower energies. However, an accu- 
rate inclusion of the such secondary emission requires de- 
tailed Monte-Carlo simulations of the EBL-induced cas- 



cades and is outside the scope of this paper. We plan 
to return to this problem in a future publication (T. 
M. Venters 2009, in preparation). Furthermore, while 
the Fermi Gamma-ray Space Telescope has begun tak- 
ing data and has already observed many more blazars 
than EGRET did, a complete catalog of blazars has not 
yet been released by the Fermi Collaboration. Thus, the 
determination of the blazar GLF(s) from Fermi data is, 
as yet, premature. With this in mind, rather than mak- 
ing a solid prediction for the anticipated Fermi results, 
we have demonstrated the sensitivity of the shape of the 
spectrum suppression due to EBL absorption to various 
model inputs, thus making explicit the information con- 
tent of such a feature. 
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spectral indices of the HBLs could possibly indicate the greater im- 
portance of external Compton emission relative to the synchrotron 
self-Compton emission that could dominate the high-energy emis- 



sion in blazars with softer spectral indices, such as the FSRQs. 



REFERENCES 



Abdo, A. et al. (Fermi LAT Collaboration) 2009, ApJ, 700, 597 
Ando, S., Komatsu, E., Narumoto, T., & Totani, T. 2007, 

Phys. Rev. D, 75, 063519 
Bottcher, M. 2007, Ap&SS, 309, 95 
Bruzual, A. G. & Chariot, S. 1993, ApJ, 405, 538 
Chiang, J., Fichtel, C. E., von Montigny, O, Nolan, P. L., & 

Petrosian, V. 1995, ApJ, 452, 156 
Chiang, J., & Mukherjee, R. 1998, Astrophys. J., 496, 752 
Chen, A., Reyes, L. C, & Ritz, S. 2004, ApJ, 608, 686 
Coppi, P. & Aharonian, F. A. 1997, ApJ, 487, L9 
Dermer, C. D. 2007, ApJ, 659, 958 

Franceschini, A., Rodighiero, C, & Vaccari, M. 2008, A&A, 487, 
837 

Giommi, P., Colafrancesco, S., Cavazzuti, E., Perri, M., & Pittori, 

C. 2006, A&A, 445, 843 
Gilmore, R. C, Madau, P., Primack, J. R., Somerville, R. S., & 

Haardt, F. 2009, arXiv:0905.1144 
Hartmann, R. C, et al. 1999, ApJS, 123, 79 
Hauser, M. G. & Dwek, E. 2001, ARA&A, 39, 249 
Kazanas, D., & Perlman, E. 1997, ApJ, 476, 7 

Kneiske, T. M., Bretz, T., Mannheim, K., & Hartmann, D. H. 2004, 
A&A, 413, 807 

Kneiske, T. M., Mannheim, K. 2005, Proceedings of the 29th 

International Cosmic Ray Conference, 4, 1 
Miicke, A., & Pohl, M. 2000, MNRAS, 312, 177 
Mukherjee, R., & Chiang, J. 1999, Astroparticle Physics, 11, 213 



Narumoto, T., & Totani, T. 2006, ApJ, 643, 81 

Padovani, P., Ghisellini, G., Fabian, A. C, & Celotti, A. 1993, 

MNRAS, 260, L21 
Pavlidou, V., Siegal-Gaskins, J. M, Fields, B. D., Olinto, A. V., & 

Brown, C. 2008, ApJ, 677, 27 
Pavlidou, V. & Venters, T. M. 2008, ApJ, 673, 114 
P rimack, J. R., G ilmore, R. C, & Somerville, R. S. 2008, 

larXiv:0811.32^0l 
Salamon, M. H., & Stecker, F. W. 1994, ApJ, 430, L21 
Salamon, M. H., & Stecker, F. W. 1998, ApJ, 493, 547 
Siegal-Gaskins, J. M., & Pavlidou, V. 2009, Phys. Rev. Lett., 102, 

241301 

Sikora, M., Blazejowski, M., Moderski, R., & Madejski, G. M. 2002, 
ApJ, 577, 78 

Stecker, F. W., Salamon, M. H., & Malkan, M. A. 1993, ApJ, 410, 
L71 

Stecker, F. W., & Salamon, M. H. 1996, ApJ, 464, 600 

Stecker, F. W., Malkan, M. A., & Scully, S. T. 2006, ApJ, 648, 774 

— . 2007, ApJS, 658, 1392 

Strong, A. W., Moskalenko, I. V., & Reimer, O. 2004, ApJ, 613, 
956 

Venters, T. M. & Pavlidou, V. 2007, ApJ, 666, 128 
Venters, T. M. 2009, in preparation 



APPENDIX 

A. SAMPLE BIAS CORRECTION M(a) 

The ISID p(a) obtained by VP07 using EGRET data is measured using a roughly flux-limited data set (to the extent 
that we can postulate that EGRET resolved all blazars of integrated gamma-ray flux greater than some value -F 7 , m i n 
and none with smaller F 7 ), so that 

1 f°° d 2 N 
p(a) = — / -TFnr^ (Al) 



f =f min dF^da 
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where AT to t is the total number of objects in the sample. However, this is not the SID pz,(a) that enters Eq. [8j The 
latter is defined by Eq. (fil) and is the distribution of spectral indices for blazars in a luminosity interval between L 1 
and L 7 + dL 1 . A flux- limited sample will be biased toward harder spectral indices than a fixed gamma-ray luminosity 
interval, because not all blazars with the same L 7 have the same flux: harder blazars have higher fluxes in the high- 
energy band and are more easy to detect. A relation can be derived between p(a) and pl(oi), starting from a relation 
between d 2 N/dF^da and d 3 'N / l dL 1 dV cora da: 



d 2 N 
dF 1 da 



dz 



d 3 N 



z=0 



dLjdV c( 



d(L 7 , V com ,a) 



d(F 7 , z, a) 



(A2) 



L 7 is proportional to F-y multiplied by a function of z and a (see Eq. Bll. In addition, in writing Eq. ([!]), we have 
assumed that a is independent of L 7 and z. Thus, we obtain 



d(L 7 , V" com ,a:) 
d(Fry, z, a) 



<9L 7 dV r 



dF^ 



com 

dz 



i 
F, 



dV r 



com 

dz 



since 



L 7 = AndKa - 1)(1 + z) a ~ 2 E f F 1 . 



Eq. (A3 1 combined with Eqs. (A2| and Q, gives 

d 2 N 
dF-d^ =PL{a) 

Substituting into Eq. ( |Al| we obtain 

1 



1 f°° 

— / dzL 1 p 1 (L 1 

?1 Jz=0 



dV C( 



dz 



p(a) 



Pl{o) 



dF^ 



dV r< 



7 J z=0 



dz 



; p L (a)M(a) 



where the last equality gives the definition of the sample correction bias M(a); the normalization, N tl 
by requiring that pl(oi) integrates to 1. 



(A3) 
(A4) 

(A5) 

(A6) 
is obtained 



B. DERIVATION OF THE BLAZAR CONTRIBUTION TO THE EXTRAGALACTIC GAMMA-RAY BACKGROUND 

The isotropic gamma-ray luminosity of a blazar at some fiducial rest-frame energy, Ef (the energy emitted in photons 
of energy Ef per unit time, assuming that the blazar emits isotropically) , is related to its integrated photon flux, F^ 
(the number of photons emitted in energies above observer frame energy Ef per unit time per unit area), through 



47rd£(a - 1)(1 + z) a - 2 E f F 1 , 



(Bl) 

where d^ is the luminosity distance, and we have assumed that the blazar has a single-power-law energy spectrum 
{dN 1 /dE 1 oc E~ a ). In turn, the differential single-blazar flux, -Fph.i^o) (number of photons per unit energy per unit 
time per unit area emitted at observer- frame energy, Eq), is related to F~ through 



F 7 = 



^ph,i(^o) dE . 



Having assumed a power-law spectrum, F 1 becomes 



F^ / F ph>1 (E f ) 



E 



E 



dE 



EfF phA (E f ) 



1 



Substituting the above into the equation for L 1 and solving for F p h,i(-E/), we get 



F phA (E f ) 

Thus, neglecting absorption, Fp^^^Eg) is given by 

E p h,i(Eo) — -F p h,i(^/) 



4ird 2 L E 2 f 



Eo 
Ef 



\itd\E) 



Eo 
Ef 



Including the effects of absorption, we arrive at Eq. ([6| of Section 2 

L. 



E p h,i(E ) — - 
The intensity of emission is defined as 

Ie(E) 



t{1 + zY 



l L^f 



Eo 
Ef 



dtdAdEdn 



_d_ 

dn 



dNF, 



(B2) 



(B3) 



(B4) 



(B5) 



(B6) 



(B7) 
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where N y is the number of photons, N is the number of objects, and -Fph.i is the flux from a single contributing object. 
Making the dependencies explicit, the intensity can be expressed as 

d f d 3 N 
Wa) = dnj Wfl., *, L„a) dL ^ dVcomda dL 7 dV com da. (B8) 

The differential number of objects can be expressed in terms of the GLF and the ISID: 

d 3 N 

= p 1 {L 1 ,z)p L {a). (B9) 



dL^dV cora da 



The comoving volume element is given by (assuming ACDM cosmology) 

d 2 V r 

11 'com 



dzd n =^^ + ^(1 + (bio) 

where D — + z) is the distance measure. Substituting for F p h.i, dV com , and the differential number of objects, 

we finally arrive at the sought for expression for the emission without attenuation: 

I E (E ) = ^4^2/ dap L (a) (J) Jdz (1 + z)~ a [Q A + Q m (l + zf]- 1 ' 2 J dL^p^ (Bll) 

The above expression for the collective unresolved blazar emission is equivalent to Equation (10) of SS96, corrected 
for cosmology and assuming that blazars consist of a single spectral population. 
Including the attenuation factor, exp [—t(E , z)], in the expression for the single-blazar flux, we obtain: 

Ie(E ) = Jf o ^E}J dap L (a) (J^j J dze-< E ^\l + z)~ a [n A + Q m (l + z) 3 }^ 2 J dL 7 L 7 p 7 . (B12) 



